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A REVIEW  OF  THE  STATUS  OF 


FIRE  PREFORMANCE  PREDICTIVE  METHODS 
by 

Stanley  P.  Rodak;  Physicist,  Heat 


The  problems  of  fire-resistance  have  attracted  the  attention  of 
various  workers.  This  project  is  concerned  with  the  theoretical  aspects 
of  fire  resistance.  From  a knowledge  of  the  thermal  propertities  of  a 
building  material  shaped  into  a particular  geometry,  we  wish  to  be  able 
to  predict  the  thermal  fire  rating**of  scaled  constructions  of  the  same 
material.  Ideally,  we  would  like  the  calculations  techniques  to  include 
variable  thermal  properties  and  endothermic  and  exothermic  processes. 


In  order  to  circumvent  certain  numerical  difficulties  in  using  vari- 
able thermal  properties,  it  was  felt  that  sufficient  information  would 
be  had  if  the  calculations  are  made  in  a manner  that  will  allow  the  infor- 
mation to  be  displayed  in  the  following  fashion.  Giedt  (1)  has  suggested 
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a mean  value  for  thermal  conductivity  when  it  varies  with  temperature. 


Endothermic  and  Exothermic  processes  are  to  be  disregarded  until 
the  numerical  techniques  are  sufficiently  developed. 


The  mathematical  equations  to  be  used  are  as  follows: 


parabolic  equation  of  heat  flow* 

-Q||)  T = 0 

(1) 

boundary  condition 

, St 
^Sn 

= -hT 

(2) 

(n  is  the  normal  surface  component) 


* Nomenclature  is  at  end  of  text. 

**  The  thermal  fire  rating  will  be  taken  as  the  time  before  the  tem- 
perature in  a region  in  the  specimen  under  consideration  reaches 
a value  predetermined  to  be  unsafe  for  fire-resistance  purposes. 
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The  mathematical  model  which  we  will  work  with  is  a homogeneous 
solid  bounded  by  two  parallel  planes.  The  problem  will  be  two-dimensional. 

Over  one  surface,  the  exposed  surface,  we  define  a temperature  rise 
by 


T(o,y,t)  = f^  (t)  (3) 

where  f^  (t)  is  the  ASTM  time-temperature  curveo 

On  the  other  surface,  the  unexposed  surface,  the  heat  flux  is  that 
for  a horizontal  surface  (2,3) 

h = h + h (4) 

n c 

h = (.174e/(T  -T  ))  • ((T  /lOO)^  - (T  /lOO)^)  (5) 

XI  so  s o 

1/3 

h = .275  T ' (6) 

c s 


For  a high,  vertical  surface,  h c would  still  be  proportional  to  T 
but  the  coefficient  changed  (3,4).  The  mechanism  of  transpiration 
cooling  was  not  considered  (5). 


The  numerical  technique  we  choose  should  be  accurate  enough  to  calcu- 
late temperatures  for  our  test  example  (6,  pg  126). 

TEST  CASE: 


The  region  o<x<i  has  zero  initial  temperature.  The 
end  x=o  is  kept  at  temperature  To  for  t>o.  At  x-X  there 
is  radiation  into  a medium  at  zero  temperature;  then, 
temperature  for  o<x<i.  are  given  by 


T(x,£)  = T 


l+i'h(i-x/i') 

I + ih 


Li 


2(Pn  + (^h)^  sin  (3n 

\n  (Xh+(lh)"+rn) 


-Pn 


(7) 


where  Pr  are  the  positive  roots 

The  exponential  term  causes  the 

Three  finite-difference  methods 
ferential  equations  have  been  studied 


of  Pn  = 0, 

series  to  converge  rapidly, 
of  approximating  the  partial 


dif- 
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The  first  studied  was  the  forward-difference  explicit  algorithm  (8); 


pcAx*" 


(T.  . 

^,3 


- .)  = kAt  (T^  + - 4T^  .) 

i^j  1-1 J 1+1;  J i;J-l  ^^J+^ 


(8) 


Here,  AY  ..  AX 

For  a slab  NAX  thick  and  MAY  tall, 

Tj  = f (kAc)  J = 1,  2,  ...,  m (9) 

j 

t lx 

is  the  temperature  on  the  exposed  surface  at  the  k interval  of  time*. 

The  requirement  that  At /Ax“  be  chosen  such  that 

kAt  ^ I 

pcAx^  4 (10) 

in  order  to  have  stability,  that  is,  the  finite-difference  algorithm  is 
immune  to  the  accumulation  of  round-off  errors  (7),  caused  a survey  of 
other  methods  (metal  rods  were  to  be  imbedded  in  the  test  example  interior). 

The  backward  difference  implicit  algorithm  (8)  is  defined  by 


pcAx^  (T^'^J 

J 


J 


.)  = kAt  (T 


k+1 

i-1 


,3 


+ T 


k+1 
i+1,  j 


+ T 


k+1 


i,j-l  i,j+l 


k+1 

- 4t7!)  (11) 


This  method  is  stable  for  all  values  of  At/Ax'".  Iterative  and  matrix 
reduction  techniques  are  available  to  solve  the  resulting  unknown  equa- 
tions (9,10,11,12).  For  a solid  of  N interior  points  in  each  direction, 
equations  result  for  the  two-dimensional  problem. 

The  backv/ard-dif f erence  implicit  algorithm  was  applied  to  a parti- 
cular problem:  finite  thick  slab  exposed  on  one  surface  to  temperature 


* 


Ic  Ic 

To  meet  the  boundary  condition  (2),  we  replace  (T . ^ . - T,  , ) by 

Xj  - 

k i 

T^  j)  by  1/2  in  equation  8.  The  unexposed  surface  is  at  i = + 1; 


i 1 is  the  exposed  surface,  we  make  a similar  replacements  in  (11) 
and  (12)  to  meet  boundary  conditions  of  (2). 
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unexposed  surface  has  constant  "cooling”  coefficient  ho,  metal 
rods  were  placed  1"  from  exposed  face.  The  equations  were  solved  for 
each  successive  time  increment  (At  = 1 min..  Ax  = Ay  = .2")  by  Gauss 
Siedel  Iteration*.  Figs.  1 & 2 are  graphs  of  the  computer  solved 
problem . 

Fig,  2 contains  two  different  profiles  of  the  distance  from  exposed 
slab  surface  versus  temperature.  One  profile  is  taken  along  a line  (AA ' ) 
normal  to  the  slab  surface  and  passing  through  one  of  the  steel  rods. 

The  other  is  along  a line  (BB*),  parallel  to  AA',  but  about  3"  from  the 
steel  rod  (the  steel  rods  imbedded  in  the  concrete  slab  were  on  6"  centers) 
Temperatures  at  30  min.,  50  min.,  and  210  min.,  intervals  along  the  line 
AA ' from  0 to  .8"  have  negative  curvature  (13).  This  is  not  the  correct 
slope  ( 14) . 

The  Crank -Nicolson  algorithm  (15)  is  reported  to  be  unconditionally 
stable  (16).  It  is  obviously  based  on  sounder  physical  arguments  than 
the  two  previous  methods.  The  algorithm  is 


. 2 /„k+l 
pc Ax  (T , , 

J 


t;  ) = (kAt/2) 
J 


( k+1  k+1  k+1  k+1 

^'i-1,  j ^ \+l,  j ^ j-1  ^ j+1  ^ 


T^  + T*f  . + T^  , T + T^f  - 4T^"^^  - 4t'!^  .) 

1-1,  j 1+1,  j 1,  j-1  i,j+l  i,j  i.j 


(12) 


The  same  test  case  (metal  rods  imbedded  in  a slab  as  used  for  the  backward 
difference  implicit  algorithm  was  used  to  test  this  algorithm.  Graphs  of 
the  computer  solved  problem  are  displayed  in  Fig.  3 & 4 . The  distance- 
temperature  curves  in  Fig.  4 have  the  proper  curvature. 


As  a precursory  check  of  the  Crank-Nicolson  algorithm,  it  was  used 
to  calculate  surface  temperatures  of  the  test  example  (homogeneous  solid 
bounded  by  parallel  planes).  The  resulting  surface  temperature  calcula- 
tions, as  well  as  those  predicted  b}’'  theory,  are  displayed  in  Fig.  5. 

The  surface  temperatures  differ  by  307o  at  100  minutes. 

In  order  to  make  the  Crank-Nicolson  difference  approximation  more 
closely  approximate  the  differencials,  equations  were  derived  that  allow 
the  Ax  spacing  to  vary  in  an  unequal  manner  through  the  slab;  the  spacing 
Ay  is  equal.  This  would  allow  a fine  grid  near  the  unexposed  surface, 
where  the  temperatures  vary  more  slowly  with  time.  Note  2 discusses  how 
a grid  spacing  was  chosen  for  .17,  "accuracy." 


Another  problem  was  encountered,  for  the  above  example  used  as  a 

check  in  the  variable  Ax  computer  program,  interior  temperatures  were 

higher  than  surface  temperatures  , k i nr.  i i • i o 

(T  = 100  F,  k = I,  . . . ; J = I,  2,  . 

J 

* The  equation  showed  diagnol-dominance,  hence  the  choice  of  Gauss- 
Siedel  Iteration  techniques.  See  Note  1 for  termination  procedures. 
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1 2 

The  aleorithm  values  had  definitely  conv®rged(T,  . < T„  . ~ 139.9  F 

1,J  2,  j 

after  60  iterations).  The  conditions  under  which  this  algorithm  will 
converge  to  values  greater  than  the  exposure  surface  temperatures  (and 
hence  apparently  violate  the  first  law  of  thermodynamics)  are  discussed 


2 2 

in  Note  3.  This  discussion  is  for  T ^ i where  the 

2 

tion  quess”  of  T.  . = 0 i = 2,  Xoj  j = 1,  2 ...^ 

j 

2 2 

found  that  similar  unstable  conditions  (i.e.  T.  . < T_  ) 

J • J 

2 2 

if  the  initial  iteration  quess  was  T.  , = T.  .or  even  if 


"initial  itera- 

Yq.  We  also 

occurred  even 

T^  . 

j were 


assigned  their  exact  values  as  predicted  from  theory  (the  convergence 


of  T 


2 

2.j 


was 


still  140  °F). 


The  present  status  of  the  project  is  to  circumvent  this  problem.  It 
may  be  that  the  constraints  involved  in  keeping  the  Crank -Nicolson  algorithm 
from  violating  the  first  law^  a At  several  orders  of  magnitude  smaller  will 
have  to  be  used. 
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NOMENCLATURE* 
surface  film  coefficient 

convection  heat-transfer  film  coefficient 

radiation  heat  transfer  film  coefficient 

conductivity 

thickness 

temperature 

unexposed  surface  temperature 
exposed  surface  temperature 

finite-difference  notation  for  temperature  at  time 

increment,  x coordinate  of  iAx  = Ax  and  y coordinate 
jAy  ..  jAx. 
t ime 

space  coordinates 
dif fersivity 
surface  emissivity 
density 

finite  division  of  x or  y coordinate 
finite  division  of  time 


* 


Units  are  English 
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Note  1 


A DISCUSSION  OF  WHEN  TO  TERMINATE  CALCULATIONS 
k+1 

OF  T . BY  GAUSS-SIEDEL  ITERATION* 

J 


Consider  the  equations 


(4  + ^7^)  . + T^^J  . + 1 + T^"*"!  . + . 

kAt  1,.]  1-1,1  i+l,J  i,J+l  J-1  i,j 


(1) 


+ Ax  ^ + 2)  T^"-^  = ^ (T^-"!  . + . 

kAt  k i,j  1-1,  J 2 i,J-l  i,J+l  i,J 


(2) 


where 


a°i,j  = T^^  . + t;  . + t.  + T‘f  . , - 4T, 

' 1-1, J 1+1, J ■ 


^ + T^  - 4T^ 

iJ+1  iJ-1  i,j 


a"  i,  j = TV  . . + ^ (T^^  . , + T;'  ) + ( 

' 1-1,  J 2 i,j-l  1,  J+1 


pc  . , h 


-Ax  f -2)T?  . 
kAt  k 1,  j 


Equation  1 is  the  Crank -Nicolson  (CN)  algorithm  for  finding  interior 
temperatures  of  a solid  bounded  by  two  infinite  parallel  planes.  Equa- 
tion 2 is  the  CN  algorithm  for  the  unexposed  surface  temperatures.  (The 
solid  is  NAx  thick  and  MAy  tall  (Ax  = Ay)).  Since  the  equations  are 
diagonally  dominate,  Gauss-Siedel  Iteration  techniques  are  used  to  solve 

the  N'M  unknown  „k+l  ..  m.i  , too  • , • 

T^  ^ (i  = 2,3,  ...,  N+1;  j = 1,2, 3 M;  i = 1 is 

the  unexposed  surface). 

k k+1 

All  T.  , and  T.  . are  known. 
i,J 

k+1 

As  an  initial  start  to  the  iteration,  all  T.  , are  set  equal  to 

' 1.  J 

2 2 

zero.  Equation  1 is  used  to  calculate  T„  ^ . This  new  value  of  T„ 

z,  1 2,1 

2 

is  used  now  in  the  calculation  of  T„  The  iteration  process  of 

k+1 

replacing  the  old  values  of  T^  . by  the  newly  calculated  values  of 
k+1 

T.  . continues  methodically  by  working  through  the  j = 1,2,  ...  M for 

J 


* The  method  of  determining  the  termination  of  the  calculations 
is  mainly  from  (21) . See  also  (9) . 
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k+1 

each  i index  until  a new  set  of  T , , have  been  iterated. 

J 

Let 

n k+l  lc+1 

AT , , = (T  (nth  iteration)  - T.  . (n-  1 iteration) 

(In  the  following  discussion^  it  is  convenient  to  consider  only  values 
of  n > 3 . ) 


For 


our  diagonally  dominate  matrix.  AT?  . should  become  small 


after  a finite  number  of  iterations.  That  is^  we  assume 


lim  at:  . = 0 

J 


n-e  eo 


i = 2,3,  . . .,N 
j = 1,2,...,M 


In  practice,  this  process  must  be  terminated  after  a finite 

number  of  iterations  (Iq)*  This  is  called  truncating.  Let  S.  . 

J 

be  the  '"exact"  temperature  value  at  [(i+l)Ax,  JAy].  The  radius  of 

k+1 

convergence  is  given  by  (T . . in  this  formula  assumes  the  value  of  the 

io  iteratioii) : 


r.k+1 

T.  . - S.  . 
J iJ 


< 


1> 

1-L 


AT^  . 

^,2 


where* 


L = sup 
n=3 


sup 
i = 

j = 1,  2. . o,M 


i = 2,3.  . .,N 


j at’^'^^: 

1 — 

I AT^ 

' J 


That  is,  for  the  (i+1)  iteration,  we  find  the  maximum  value  \t  of  the 
-T+1  i 


ratio 


AT 


u. 


AtJ  . 1 

J 


A._  = sup 


i = 2,3, 

j - 1.2, 


M 


A T + 1 

AT.  . 


ATT  . 
J 


reference  13,  pg  9 
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From  the  set  of  (3  < T < <»),  we  again  find  the  maximum  value  L, 

OO 

L = sup  \ 

T..3 


In  practice^  an  approximate  value  of  L is  obtained  from  the  set 
of  \ of  the  first  ten  iterations., 

T 


Note  2 


GRID  SIZE  DETERMINATION 


The  following  houristic  method  of  arriving  at  an  equation  to 
determine  the  proper  grid  size  to  use  in  one-dimensional  finite  dif- 
ference equations  was  described  by  Dr.  Arms  (21).  We  can  extend  the 
results  to  two-dimensional  problems.  Let  f.  ((i+1)  Ax)  be  the  itera- 
tived  temperature  value  for  grid  size  Ax.  Call  f((i+l)  Ax)  the  "exact" 
value.  Then^  in  an  approximate  manner^ 

£Ax  ((i+l)^x)  = f((i+l)  Ax)  + Ax^.A 


A can  be  thought  of  as  an  operator.  (The  series  expansion  was  terminated 
on  the  first  term).  Taking  a grid  size  of  Ax/2 


'Ax/ 2 


= f + 


Ax^ 


or 


f . - f Ax/2 

Ax 


A + 2nd  order  terms. 


By  doing  the  same  problem  for  two  different  grid  sizes^  A can  be 
determined  (we  neglect  the  2^^^  order  terms),  and  hence  the  proper 
variable  grid  spacing  for  the  desired  degree  of  accuracy. 


- s*.- . 


Note  3 


A DISCUSSION  OF  THE  CONSTRAINTS 
IMPOSED  ON  CRANK-NICOLSON  FINITE -DIFFERENCE 
ALGORITHM  BY  GAUSS -SIEDEL  ITERATION 


Consider  a solid  (i  thick)  bounded  by  two  parallel  planes.  On 
one  surface^  there  is  a sudden  temperature  rise  Ts  for  t>0.  The  ini- 
tial temperature  of  the  solid  is  zero.  Using  Crank-Nicolson  difference 
equations^  we  wish  to  find  the  interior  temperatures  at  time  At.  Gauss- 
Siedel  Iteration  will  be  used  to  solve  for  T^  • . As  a first  quess  set 

all  T?  . equal  to  zero.  The  equations  for  a irariable  Ax  and  Ay  are: 

J 


i=l 


exposed 

surface 


for  t>o) 


inside  solid: 


(9sAll.  + ^ + AiL + 1) 

^kAt  Ax^  (Axj^+Ax^)  Ax2(Ax^+Ax2>  i,  j 


^ 

Ax^(Ax^  + AX2) 


k+1  Av^ k+1  1 

i-l,j  Ax2(Axj^  + Ax^)  i+1,  j 2 


J-1 


+ 


T 


k+1 

i,  j+1 


) 


(1) 


_ Ay^ k Ay^ k+1  1 .k+1 

^i,j“  Axj^CAy^  + Ax2>  Ax2(Ax^  + Ax2>  i+1,  j 2 ^ i,  j-1 

) 4f£gAyl ^ k 

i^  j+2^  \kAt  Axj^(Ax^  + Ax^)  Ax2(Ax^  + Ax^)  / ^i 
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at  surface  'k  = i 


i+l,j 


Axj  k ^ kA:  ^ ^i,J 


/AZ  N2  ™k+l 
'^Ax/ 


, 1 /„k+l  , „k+l  . , 

i-lj  2 ^i,j+P  ■*■  \ j 


^kAt 


- ( 


Ax, 


Ay  h 
Ax  k ' 


1)  T 


let  Ax-i  = Ax,  Ay  i = 2^  • * 
■1  1 


m. 


Equation  (1)  then  becomes  (noting  for  our  case  = 0); 


i,  j ^ "i+1  ^ ^ ^^i,j-l  i,j+r 


Where  0 = ( 


2pcAy'‘ 
k At 


-1 


+ 2(gr  +2)  ^ < 1 


Let  us  introduce  some  symmetry  into  the  boundry  condition: 


_k  k 

^i,  1 - ^i,3 


Then  the  first  iteration  is: 
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T?  . = gu)T  (1  + P)  = 

Zm  y ^ 8 ^ y J- 

T2  4 = ^“^8  ^ 


2 Q ^ /I  I Q . q2  . , Qin-2.  Q „ 

T»  = P(A)T  (1  + P + P + .00  + P ) = PodT  (t — Q ) 

2,  m s s i-p 

2 2 2 
t:^  „ = pVt 

3,  2 8 

„ = pVt  (1  + 2 P) 

3, 3 s 

„ = pVt  (1  + 2p  + 3P^) 

J j s 


T?  = P^OD^  n + 2P  + 3P^  + 
3,m  s 


+ (m-1) 


l-p”~^-(m-l)p"^ 

- (1  - P)^ 


etc , 


second  iteration: 


- = P [u)T  + P^U)\  + 2Puyr  (1  + p)1 
2,  2 s s s 


= PtuT  (+  2P  + 2P^  + P^oo^)  < 1 

s 


For  the  various  heat-conclusion  problems  under  study^  this  type 
of  analysis  was  used  to  understand  the  constraints  imposed  on  the  iteration 
techniques.  It  was  found  that  when  the  constraints  were  violated,  there 

was  convergence  of  the  (all  equations  had  diagonal  dominance),  but 

J 

k+1  k 

the  j were  larger  than  for  several  Ax^  from  the  exposed  surface  ^). 
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